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Abstract 

We introduce an integrable, four-well ring model for bosons where the 
tunneling couplings between nearest-neighbour wells are not restricted to 
be equal. We show how the model may be derived through the Quantum 
Inverse Scattering Method from a solution of the Yang-Baxter equation, 
and in turn solved by algebraic Bethe Ansatz means. The model admits 
multiple pseudovaccum states. Numerical evidence is provided to indi¬ 
cate that all pseudovacua are required to obtain a complete set of Bethe 
eigenstates. The model has the notable property that there is a class of 
eigenstates which admit a simple, closed-form energy expression. 
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1 Introduction 


Systems of Bose-Einstein Condensates (BECs) continue to offer rich opportu¬ 
nities for exploring intrinsically quantum characteristics, such as superposition 
and entanglement, at the mesoscopic/macroscopic interface. Models based on 
few bosonic modes have been widely studied. Initially, a two-mode model of 
quantum tunneling between two wells, specifically the two-site Bose-Hubbard 
model, was sufficient to predict self-trapping phenomena hhh which was subse¬ 
quently confirmed experimentally [5]. Ensuing studies were extended to three¬ 
mode models to investigate, for example, dynamics (BU7| . entanglement, |S], dis¬ 
sipation [5] , quantum phase transitions urn. and condensate fragmentation mi- 

The dynamics of a four-well Bose-Hubbard model was studied in [T2] as a 
means to achieve mass transport and persistent currents around a loop. The 
same model, but with two significantly different tunneling rates, has been in¬ 
vestigated in relation to a mesoscopic quantum system in thermal contact m- 
The quantum dynamics for a range of different initial conditions, in terms of 
the particle number distribution among the wells and quantum statistics, was 
presented in |14l . In |15j it is claimed that an appropriate control of short-range 
and dipolar interaction may lead to the dynamical creation of mesoscopic quan¬ 
tum superposition, which could be employed in the design of Heisenberg-limited 
atom interferometers through a four-well system. Also, it has been proposed to 
implement a two-well model with two levels in each well (yielding a four-mode 
model) as a means to experimentally measure EPR entanglement |16j . 

The two-mode Bose-Hubbard model is an example of an integrable system 
which admits an exact Bethe Ansatz solution ffZHISI, the existence of which 
allows analytic computation of physical quantities. For example, a Bethe Ansatz 
solution of the model was used to calculate the ground-state one-body density 
matrix in the attractive regime [20] . This calculation shows the existence of a 
quantum phase transition point of the model which is characterized as separating 
condensate fragmentation from an unfragmented phase. It is desirable from this 
perspective to also identify integrable cases of similar models with higher number 
of modes. This is the objective of our work here for the case of a four-well model. 

In Sect. 2 we define the model and set up the notational conventions. Sect. 
3 is devoted to deriving the model via the Quantum Inverse Scattering Method, 
and deriving an exact solution through the algebraic Bethe Ansatz mu- Here 
we encounter some unusual properties for this model. One is that although we 
can identify four conserved operators for the system, only two of them are ob¬ 
tained through the transfer matrix produced by the Quantum Inverse Scatter¬ 
ing Method. Another curiosity is that there are multiple pseudovacuum states 
available for the purposes of implementing the algebraic Bethe Ansatz. A con¬ 
sequence of this is that there is a class of eigenstates which admit a simple, 
closed-form energy expression. Numerical studies undertaken in Sect. 4 indi¬ 
cate that all pseudovacua are required in order to obtain a complete set of Bethe 
eigenstates. Concluding comments are given in Sect. 5. 
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2 Hamiltonian of a four-well ring with anisotropic 
tunneling 

The four-well model we present has the Hamiltonian 

H = U(Ni + N 3 -N 2 - N 4 ) 2 + ii(Ni + N 3 - N 2 - N 4 ) 

+ ti 2 (aia\ + a\a 2 ) + *14(0104 + a\a 4 ) (1) 

+ *23(0303 + 0302) + *34(0304 + aga 4 ). 

The operators a] , a t , i = 1,... ,4 denote the single-particle creation and annihi¬ 
lation operators acting in the four wells, with iVj = a\ai providing the number 
operator for bosons in well i. The quantities which determine the tunneling 
couplings are not independent but take a factorised form tij = —where 
the set {ft, aj, j = 1,2,3,4} are arbitrary real variables. This is equivalent to 
the constraint 


*12*34 — *23*14 

but still admits sufficient freedom to investigate a range of anisotropic tunnel¬ 
ing regimes. The coupling parameters [/, which controls on-site and inter-well 
interactions between bosons, and /i, which is an external potential parameter, 
are also arbitrary real variables. 

The total number operator N = 1 *Vj is a conserved operator commuting 

with the Hamiltonian. For each fixed value of N the dimension of the associated 
Hilbert space Vn is 

dim (Vat) = + 3 )(N + 2 )(N + 1). 

Later will show that there are two additional constants of the motion, which 
leads to the conclusion that the Hamiltonian is integrable. In Fig. [T]we show a 
schematic representation for the model. 



Figure 1: Schematic representation of the model m showing the tunneling 
couplings among the wells. 


Physically, the Hamiltonian m models anisotropic Josephson tunneling for 
a Bose-Einstein condensate confined in a four-well ring potential. This Hamil¬ 
tonian differs from the usual Bose-Hubbard model mm by the inclusion of 
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interactions terms of the form NiNj. However these quadratic interactions are 
found in other models, e.g. mm- The inclusion of the external potential term /z, 
which has an analogue in the case of the two-mode Bose-Hubbard model pToTF27j] , 
is compatible with the integrability. 

There are rich mathematical structures associated to this model, including 
an exact Bethe Ansatz solution which will be derived below. 


3 Quantum Inverse Scattering Method 


Our first objective is to establish that the Hamiltonian CD can be obtained 
through the Quantum Inverse Scattering Method. We follow the general ap¬ 
proach of [2BJ, which itself extends ideas presented in m- We begin with the 
standard sit(2)-invariant I?-matrix, depending on the spectral parameter u: 


/ 1 0 0 0 \ 

0 b(u) c(u ) 0 

0 c(u) b(u) 0 

\ 0 0 0 1 / 


( 2 ) 


with b(u) = u/(u + rf) and c(u ) = 77 /( 1 / + rf). Above, r) is a free real parameter. 
It is well-known that R(u) satisfies the Yang-Baxter equation P51I25] 


Ri 2 (u - v)R 13 (u)R 23 (v) = R 23 {v)Ri 3 {u)Ri 2 (u - v). ( 3 ) 


Here Rjk(u) denotes the matrix acting non-trivially on the j -th and fc-th spaces 
and as the identity on the remaining space. We start with the Lax operator m 


L {i ’ j) (u) = 


ul + TjNi'j 

At. 




V 1 ( a i + °^)i. 


(4) 


where I denotes the identity operator and 

A Mi -f- Mj , A,;j ck'/fij (_X'j dj , A if - (y t e, -f- otjCtj. 

Since 1,0b) (u) satisfies the Yang-Baxter equation 

R\ 2 (u - v)Li’ j \u)L^\v) = L 2 ’ j \v)L^ (u)R 42 (u - v), ( 5 ) 


we can define a monodromy matrix from the Yang-Baxter realization presented 
above as 

= L™(u + w)L™( U - u) = ^(“) , ( 6 ) 

where 


A(u) = ((it + u)I + r]Nij)((u - w)/ + r]N 2}4 ) + A 13 A\ a 
B( u) = ((it + ui)I + ryA r i, 3 )A 2 ,4 + + a\)A ^ 3 

C(u) = A\ 3 ((u - w)I + 77^2,4) + V "Vi + a\)A\ 4 
D(u) = A\ 3 A 24 + t] 2 (aq + a 3 )(a\ + a 2 )I. 
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It follows that the monodromy matrix satisfies the Yang-Baxter equation 

Ri 2 {u - v)Tx(u)T 2 (v) = T 2 {v)T 1 (u)R 12 (u - v). ( 7 ) 

Finally, we define the transfer matrix 

t{u) = trace(T(u)) = A(u) + D{u) = cq + c\u + c 2 u 2 , ( 8 ) 


where 


co = t(0) = —k 1 H + ^ r / 2 (a 2 + a\ )(a| + a 2 ) + - u 2 ^ I, 

ci = 4 ~ t ( u )\ u =o = rjN, 
au 

1 d 2 

C2= 2d^ r(u)l “= 0 = / 

and the following identification has been made for the coupling constants: 



/i = KU1T] tij = —KOtiOtj. 


(9) 


It follows from 0 that the transfer matrix commutes for different values of 
the spectral parameter, 

[t(u),t(v)} = 0, 

and consequently the eigenvectors of t(u) are independent of u. It is now 
straightforward to check that the Hamiltonian © is related to the transfer 
matrix t[u) © through 


H = —K 



^ 'ui 1 — u 2 — 77 2 (af + a^)(a 2 + a 2 ) — ur/N 



Therefore the energy spectrum is given by 


E = —k ( A (u) +oj 2 — u 2 — rj 2 (ay + a 2 )(a 2 + a 2 ) — ur]N — 


( 10 ) 


where A(it) denotes the eigenvalues of the transfer matrix. To find expressions 
for these eigenvalues we apply the algebraic Bethe Ansatz procedure EEJI22]- 


3.1 Pseudovacua and algebraic Bethe Ansatz 

To apply the algebraic Bethe Ansatz method we have to determine suitable 
pseudovacua. To accomplish this we define the commuting operators 

pf t t 

T[ 3 = a 3 a[ - aia£, 

pt t t 

1 2 4 — O/.4.CL2 0^2^47 
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which satisfy 

^3^1,31=0, 

P + l,3^2,4] = 0, 

[NiMrU) l ] = iFl 3 ) 1 , 

[ri,3,-B(u)] = r?ri )3 A 2 ,4, 

[ r I,3> S ( u )] = -»? r i,3^2,4, 

[r li3 ,c(«)] = 0 , 


[r^ 4 , a 2i4 ] — 0 , 

[^2,4: ^-1,3] = 0, 

[JV 2 ,4, (r+ i4 ) fc ] = fc(r+ >4 ) fc , 

[T 2A ,B(u)] =0, 
[t5 )4 ,B(u)] =0, 

[r 2 , 4) <?(«)] = f?r 2)4 A| 3 , 

p2,4iC( M )] = — ^2,4^1,3' 


Fl, 3 , <?(«)] =0, 

Let |0) denote the bosonic vacuum state defined by the property 


Oj|0) = 0 j = 1,2,3,4. 

For the purposes of the algebraic Bethe Ansatz calculation the pseudovacua are 
given by 

I0M> = ( r i, 3 ) fe ( r 2, 4 )'|0)> k + l < N (11) 

and satisfy 


A(u)| <t>k,i) = (u + u + kr]){u - w + lri)\<j)k,i), 

B{u)\<f>k,i) = 0, 

D ( u )\<f>k,i) = y~ 2 (al + a\){a\ + a\)\<j> kt i). 

Now we can define the Bethe states associated to each pseudovacuum as 

{ N—k—l 

P C(ui)\(p k ,i), if k,l = 0, l,--- ,N- 1, k + l<N 

I ill), if k + l = N. 

Using the following algebraic relations 

A(u)C(v) = - ~^-^-C(v)A{u) --— C(u)A(v), 

u — v u — v 

D(u)C(u ) = - -C(v)D(u) H--— C(u)D(v) 

u — v u — v 

which are determined from ([7|l. the algebraic Bethe Ansatz 21,^ can be ap¬ 
plied. In particular it is found that 


T(u)\l/j k ,l) = Xk,l(u)\lpk,l), (12) 

where for k + l = N the transfer matrix eigenvalues are given by 

A k,i( u ) = {u + w + krj)(u -ui + lrj) + rj~ 2 ( a l + a 2 )(a 2 + a 2 ), 
while for k + l < N the eigenvalues are given by 

N-k-l 

Xk,i(u) = (u + u + krj){u — w + lr]) TT — —— -- 

-*■ U — Vj 
3= 1 3 

N-k-l 

+ ,- 2 (a? + ^)K+oi) II "UU’’ 
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(13) 


and the parameters Vj are solutions of the Bethe Ansatz equations 

if(vi + lo + kr])(vj - lo + Irf) _ 1 Vj — Vj — r] 

(a? + «l)( a 2 + a 4 ) fj; v i - v i + V' 

These eigenvalue formulae allow for the energy spectrum to be computed through 
m- We remark that in the case when k + l = N the expression m reduces 
to the simple form 

E = —k (^klrj 2 + uir](l — k ) — ^ 

= U(l - kf + (l- k)n- (14) 

Eq. (UTl) is notable for two reasons. First it shows that there is a class of 
eigenstates which admit a simple, closed-form energy eigenvalue. The second is 
that for these eigenstates the energy eigenvalue is independent of the variables 

tij. 


3.2 Additional conserved operators 


The Hamiltonian m is a model with four bosonic modes, for which it is expected 
that integrability requires at least four independent conserved operators. How¬ 
ever, the Quantum Inverse Scattering Method applied above shows that we 
obtain only two independent conserved operators, H and N, from the transfer 
matrix. Nonetheless there exists two other independent conserved operators 
with the form 


Q 1,3 

Qza 


_pt p 

.2 1 1.3 1 1,3 


_r T p 

2 1 2,4-*- 2,4- 


These operators satisfy the commutation relations 


[Ql,3, Q2,4,] — Oi [H, Qia\ 


[H, Q 2 , 4 ] — [N, Qi,3] — [N , < 32 , 4 ] — 0 


so Qi ,3 and < 22,4 together with the Hamiltonian H and the total boson num¬ 
ber operator N provide four independent conserved operators for the model. 
Moreover, the additional conserved operators satisfy the following commutation 
relations 


[Q 1,3) C'(u)] = [<32,4, C'(u)] = 0. 

It follows that each Bethe state \ipk,i) as defined above is a simultaneous eigen¬ 
state of the additional conserved operators: 

= k\ipk,l), (15) 

<32,4|'0fc,z) = l\ (16) 


4 Numerical results for small number of bosons 

Finally we undertake a numerical analysis to demonstrate that (at least in some 
instances) utilizing all possible pseudovacua in the Bethe Ansatz procedure al¬ 
lows for a complete set of energy levels to be obtained. 
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4.1 N = 1 

By using the normalized Fock basis 

{IXj> = a ]l°) : 1 < J < 4 I 

for the sector N = 1, the Hamiltonian takes the matrix form 

(U ± /d ti2 0 tl 4 \ 

rr _ *12 U - II *23 0 , - 

0 t 23 U + IJL *34 ‘ 1 * 

\ *14 o f 34 u - n) 

Choosing the parameter values 

U = 0.125, n = -0.55, iia = f 23 = *14 = t 3 4 = -0.25, (18) 

and diagonalizing this matrix, we find that the eigenvalues are given by 

E 1 = -0.6183034374 
E 2 = -0.4249999999 
E 3 = 0.6750000000 
£ 4 = 0.8683034374 

The next step is to compute the spectrum from the Bethe Ansatz equations. 
Choosing 

= 3 = 1,2,3,4 (19) 

then © imposes the parameter values 

77=1, cj = 1.1. k = 0.5. (20) 

For the N = 1 sector there are three pseudovacua. 

• 10o,o) = |0): This pseudovacuum leads to the Bethe Ansatz equation 

r?tf-u> 2 ) = l 

which has two solutions 

V! = ±1.486606875 

Using (Hhd, this yields 

£(1.486606875) = -0.6183034375 
£(-1.486606875) = 0.8683034375 

in agreement with £i and £4 obtained by numerical diagonalization. 

• |0i, 0 ) = r{ 3 10): Here the pseudovacuum is an eigenvector of the Hamil¬ 
tonian and through m we obtain 

£ = -0.425 

which is in agreement with E 2 obtained by numerical diagonalization. 

• |0o, 1 ) = 4 10): Here the pseudovacuum is an eigenvector of the Hamil¬ 

tonian and through da we obtain 

£ = 0.675 

which is in agreement with £3 obtained by numerical diagonalization. 
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4.2 N = 2 


By using the normalized Fock basis 

{| Xi,j) = C ij a l a ]| 0 ) : 1 < * < j = 4, Cjj = 1/V2, C = 1 for * ^ }} 
for the sector N = 2, the Hamiltonian takes the matrix form 


/4U + 2 p 

H2V2 

0 

0 

0 

0 

*14^2 

0 

0 

0 

il 2 V 2 

0 

tl 2^2 

*23 

0 

0 

0 

*14 

0 

0 

0 

twV 2 

4 U - 2/i 

0 

*23\/2 

0 

0 

0 

0 

0 

0 

t23 

0 

4+ + 2p 

*12 

0 

*34 

0 

*14 

0 

0 

0 

<23\/2 

*12 

0 

*23 

0 

*34 

0 

0 

0 

0 

0 

0 

*23 V% 

4U + 2/1 

0 

0 

*34\/2 

0 

*14^2 

0 

0 

*34 

0 

0 

0 

*12 

0 

*14\/2 

0 

1 14 

0 

0 

*34 

0 

*12 

417 - 2(i 

*23 

0 

0 

0 

0 

*14 

0 

*34\/2 

0 

*23 

0 

*34\/2 

V 0 

0 

0 

0 

0 

0 

*14^2 

0 

*34 V2 

4 U - 2 11 / 


Again choosing the parameter values m to diagonalize the matrix we find that 
the eigenvalues are given by 


£i = -1.12878 
E-2 = -0.883095 
E 3 = -0.143398 

£4 = -0.600000 
£5 = 0.0000000 


£ 6 = 0.233742 
£7 = 0.2830951 
£ 8 = 1.6000000 
£9 = 1.743398 
£10 = 1.895046 


To compare these results with those obtained from the Bethe Ansatz equations, 
we again make the identifications given by uni) and For this sector we 

have six pseudovacua. 

• l^o.o) = |0): The pseudovacuum leads to the Bethe Anastz equations 



which admits the solutions 


w 2 ) = 


2 ) = 


V\ — V2 — 

V\ — V2 + TJ 
V2 ~ Vl — t] 
V 2 — Vi + 77 


{t>i = 1.128788719+ 0.4392638051*, v 2 = 1.128788719- 0.4392638051*}, 

{v! = -1.779416015, v 2 = 1.311931196}, (ui = -2.611876690, u 2 = -1.178215928} 


as well as the spurious solutions sets 

{ Vl =v 2 = ±0.4582575695} 

which must be discarded since the roots coincide (c.f. [Mj). Applying m 
we find 


£(1.128788719 + 0.4392638051*, 1.128788719 - 0.4392638051*) = -1.128788718 

£(-1.779416015,1.311931196) = 0.2337424096 
£(-2.611876690, -1.178215928) = 1.895046310 

in agreement with £ 1 , Eg and £10 obtained by numerical diagonalization. 
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|</> 1 >0 ) = r{ 3 |0): This pseudovacuum leads to the Bethe Ansatz equation 
? 7 2 (ui + lj + rj){v\ — uj) = 1. 


The solutions are 

vi = 1.386796226411321,-2.386796226411321 

which produces the energy eigenvalues 

£(1.386796226411321) = -0.14339811, 
£(-2.386796226411321)= 1.74339811 

agreeing with £3 and £9 obtained by numerical diagonalization. 

|0 o ,i) = Il 4 |0>: Associated with this pseudovacuum there is the Bethe 
Ansatz equation 

rj 2 (vi + w) (rq - ui + 77 ) = 1 (21) 

with solutions 

V! =0.666190378,-1.66619037. 

These solutions provide the energy eigenvalues 

£(0.666190378) = -0.88309518 
£(-1.666190378) = 0.2830951894 

corresponding to £2 and £7 obtained by numerical diagonalization. 

|0i, 1 ) = rJ 3 T 2 , 4 10): Here the pseudovacuum is an eigenvector of the 
Hamiltonian and through (1141) we obtain 

£ = 0 

which is in agreement with £5 obtained by numerical diagonalization. 

| 02 ,o) = (^ 1 , 3 ) | 0 ) : Again the pseudovacuum is an eigenvector of the 
Hamiltonian and through m we obtain 

£ = - 0.6 

which is in agreement with £4 obtained by numerical diagonalization. 

100 , 2 ) = (^ 2 , 4 ) |0> : Again the pseudovacuum is an eigenvector of the 
Hamiltonian and through (ThO) we obtain 

£ = 1.6 

which is in agreement with Eg obtained by numerical diagonalization. 


4.3 Completeness of the Bethe Ansatz equations 

Our numerical investigations suggest that for n = k + l < N fixed, there are 
n+1 sets of Bethe Ansatz equations m each admitting N — n + 1 non-spurious 
solutions. Moreover, when n = N there are N + 1 pseudovacuum states which 
are eigenstates of the Hamiltonian. This implies that the total number of Bethe 
states for fixed N is 

N 1 

^(?r + l)(A-n+l) = ~{N + 3)(N + 2)(JV + 1) = dim(V Ar ). 

n—0 

Thus for generic values of the coupling parameters we expect that the Bethe 
Ansatz solution of this model is complete, analogous to other exactly solved 
models where this property has been established [3D - 1M] . 


5 Conclusion 


Our main objective was to establish an exact Bethe Ansatz solution for the 
four-well Hamiltonian ©■ This was achieved by first formulating the model 
through the Quantum Inverse Scattering Method. It was found that the model 
admits multiple pseudovacuum states, and for each of these there is a set of 
Bethe Ansatz equations given by (1131) . Introducing new variables 

V = (ay + “ 3 ) 1/4 ( a 2 + a!r 1/4? ?> 

u> = {a\ + al)~ 1,4 {al + a!C 1/4 (w + - 0/ 2 )> 

Vj = {a\ + + aq)" 1/4 (^- + rj(k + l)/ 2 ) 

then m assumes the form 


fj 2 {Vi - Cj) 


N—k—l 


n 


y i - Vj -rj 
Vi - Vj + fj' 


( 22 ) 


At this point we recognise (1221) as the Bethe Ansatz equations for the two-mode 
Bose-Hubbard model UZHU. It turns out that the eigenspectrum of the four- 
well Hamiltonian © is in one-to-one correspondence with the eigenspectrum 
of a particular ensemble of two-well Hamiltonians, with particle numbers and 
coupling parameters dependent on the eigenvalues k and l of the additional con¬ 
served operators as given by m and m ■ This opens a path of investigation 
into the four-well model by utilizing results known for the two-well model. This 
approach will be developed in future work, with particular focus on character¬ 
izing the ground-state properties of the four-well model. 
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